library(readxl)
mydata <- read_xlsx("Productivity.xlsx")
mydata <- as.data.frame(mydata)
head(mydata)
##   ID Start AT1 AT2 AT3
## 1  1    85 133 136 134
## 2  2    90 123 138 135
## 3  3    96 120 135 141
## 4  4    93 129 132 134
## 5  5    84 117 138 131
## 6  6    91 128 132 141

Description:

library(rstatix)
## 
## Attaching package: 'rstatix'
## The following object is masked from 'package:stats':
## 
##     filter
library(tidyverse)
## ── Attaching core tidyverse packages ────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.2.1     ✔ readr     2.2.0
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.3     ✔ tibble    3.3.1
## ✔ lubridate 1.9.5     ✔ tidyr     1.3.2
## ✔ purrr     1.2.2
## ── Conflicts ──────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks rstatix::filter(), stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
mydata_long <- mydata %>%
  pivot_longer(cols = Start:AT3,
               names_to = "Time", 
               values_to = "Products") %>%
  as.data.frame()
  

mydata_long$Time <- factor(mydata_long$Time,
                           levels = c("Start", "AT1", "AT2", "AT3"),
                           labels = c("Start", "AT1", "AT2", "AT3"))

head(mydata_long, 10)
##    ID  Time Products
## 1   1 Start       85
## 2   1   AT1      133
## 3   1   AT2      136
## 4   1   AT3      134
## 5   2 Start       90
## 6   2   AT1      123
## 7   2   AT2      138
## 8   2   AT3      135
## 9   3 Start       96
## 10  3   AT1      120
library(ggpubr)

ggboxplot(mydata_long, 
          x = "Time", 
          y = "Products",
          xlab = "Time of measurement",
          ylab = "Products manufactured") 

library(ggplot2)
ggplot(mydata_long, aes(x = Time, y = Products)) +
  geom_violin(fill = "grey90", alpha = 0.6, trim = FALSE) +
  geom_boxplot(width = 0.1, fill = "grey60", outlier.shape = 16) +
  labs(
    title = "Distribution of Products Manufactured Across Measurement Occasions",
    x = "Measurement occasion",
    y = "Number of products manufactured per month"
  ) +
  theme_bw()

library(rstatix)
mydata_long %>%
  group_by(Time) %>%
  get_summary_stats(Products, type = "mean_sd")
## # A tibble: 4 × 5
##   Time  variable     n  mean    sd
##   <fct> <fct>    <dbl> <dbl> <dbl>
## 1 Start Products    38  89.9  4.53
## 2 AT1   Products    38 125.   4.97
## 3 AT2   Products    38 135.   5.50
## 4 AT3   Products    38 138.   4.95
library(rstatix)

mydata_long %>%
  group_by(Time) %>%
  shapiro_test(Products)
## # A tibble: 4 × 4
##   Time  variable statistic      p
##   <fct> <chr>        <dbl>  <dbl>
## 1 Start Products     0.975 0.545 
## 2 AT1   Products     0.959 0.179 
## 3 AT2   Products     0.957 0.146 
## 4 AT3   Products     0.947 0.0685
library(rstatix)

#rANOVA
ANOVA_results <- anova_test(dv = Products, #Dependent variable
                            wid = ID, #Subject identifier
                            within = Time, #Within-subject factor variable
                            data = mydata_long)

ANOVA_results #Summary of results.
## ANOVA Table (type III tests)
## 
## $ANOVA
##   Effect DFn DFd      F        p p<.05   ges
## 1   Time   3 111 760.45 8.09e-74     * 0.938
## 
## $`Mauchly's Test for Sphericity`
##   Effect     W     p p<.05
## 1   Time 0.981 0.983      
## 
## $`Sphericity Corrections`
##   Effect   GGe       DF[GG]    p[GG] p[GG]<.05   HFe       DF[HF]
## 1   Time 0.987 2.96, 109.54 7.01e-73         * 1.082 3.25, 120.11
##      p[HF] p[HF]<.05
## 1 8.09e-74         *
get_anova_table(ANOVA_results, correction = "auto")
## ANOVA Table (type III tests)
## 
##   Effect DFn DFd      F        p p<.05   ges
## 1   Time   3 111 760.45 8.09e-74     * 0.938
library(effectsize)
## 
## Attaching package: 'effectsize'
## The following objects are masked from 'package:rstatix':
## 
##     cohens_d, eta_squared, omega_squared
interpret_eta_squared(0.938, rules = "cohen1992")
## [1] "large"
## (Rules: cohen1992)
library(rstatix)

#Comparing all paires of variables
pwc <- mydata_long %>%
  pairwise_t_test(Products ~ Time, 
                  paired = TRUE,
                  p.adjust.method = "bonferroni")

pwc
## # A tibble: 6 × 10
##   .y.      group1 group2    n1    n2 statistic    df        p    p.adj
## * <chr>    <chr>  <chr>  <int> <int>     <dbl> <dbl>    <dbl>    <dbl>
## 1 Products Start  AT1       38    38    -32.4     37 9.74e-29 5.85e-28
## 2 Products Start  AT2       38    38    -39.7     37 6.28e-32 3.77e-31
## 3 Products Start  AT3       38    38    -44.7     37 8.16e-34 4.89e-33
## 4 Products AT1    AT2       38    38     -8.35    37 4.88e-10 2.93e- 9
## 5 Products AT1    AT3       38    38    -11.6     37 6.70e-14 4.02e-13
## 6 Products AT2    AT3       38    38     -2.49    37 1.74e- 2 1.05e- 1
## # ℹ 1 more variable: p.adj.signif <chr>
library(afex)
## Loading required package: lme4
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## Registered S3 method overwritten by 'lme4':
##   method           from
##   na.action.merMod car
## ************
## Welcome to afex. For support visit: http://afex.singmann.science/
## - Functions for ANOVAs: aov_car(), aov_ez(), and aov_4()
## - Methods for calculating p-values with mixed(): 'S', 'KR', 'LRT', and 'PB'
## - 'afex_aov' and 'mixed' objects can be passed to emmeans() for follow-up tests
## - Get and set global package options with: afex_options()
## - Set sum-to-zero contrasts globally: set_sum_contrasts()
## - For example analyses see: browseVignettes("afex")
## ************
## 
## Attaching package: 'afex'
## The following object is masked from 'package:lme4':
## 
##     lmer
rm_model <- aov_ez(
  id = "ID",
  dv = "Products",
  within = "Time",
  data = mydata_long
)
residuals_rm <- residuals(rm_model$lm)
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
## The following object is masked from 'package:purrr':
## 
##     some
qqPlot(residuals_rm)

## [1]  8 92